/*
 * NewtonSolver.h
 *
 *  Created on: 20 Oct 2011
 *      Author: allan
 */

#ifndef NEWTONSOLVER_H_
#define NEWTONSOLVER_H_

// Eigen includes
#include <Eigen/Core>
#include <Eigen/LU>
using namespace Eigen;

// GeoReact includes
#include <Utils/ConvergenceDiagnostics.h>

class NewtonProblem
{
public:
	virtual void Function(const VectorXd& x, VectorXd& F) = 0;
	
	virtual void Jacobian(const VectorXd& x, MatrixXd& J) = 0;
	
	virtual void FunctionJacobian(const VectorXd& x, VectorXd& F, MatrixXd& J) { Function(x, F); Jacobian(x, J); }
};

class NewtonSolver
{
public:
	enum Decomposition 
	{ FullPivotingLU, PartialPivotingLU };
	
public:
	NewtonSolver();
	
	NewtonSolver(unsigned dimension);
	
	void Initialiaze(unsigned dimension);
	
	void SetDecompositionAlgorithm(Decomposition decomposition);
	
	void SetConvergenceOptions(double tolerance, uint maxiterations);
	
	bool Iterate(NewtonProblem& problem, VectorXd& x);
	
	ConvergenceDiagnostics Solve(NewtonProblem& problem, VectorXd& x);
	
private:
	/// The residual vector	
	VectorXd F;
	
	/// The Jacobian matrix
	MatrixXd J;
	
	/// The full LU decomposition instance
	FullPivLU<MatrixXd> full_lu;
	
	/// The partial LU decomposition instance
	PartialPivLU<MatrixXd> partial_lu;
	
	/// The algorithm used for the pivoting in the lu decomposition
	Decomposition decomposition;
	
	/// The residual tolerance of the non-linear function vector
	double tolerance;
	
	/// The maximum number of iterations
	uint maxiterations;
};

#endif /* NEWTONSOLVER_H_ */
